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We present a general, systematic, and efficient method for decomposing any given exponential 
operator of bosonic mode operators, describing an arbitrary multi-mode Hamiltonian evolution, 
into a set of universal unitary gates. Although our approach is mainly oriented towards continuous- 
variable quantum computation, it may be used more generally whenever quantum states are to be 
transformed deterministically, e.g. in quantum control, discrete-variable quantum computation, or 
Hamiltonian simulation. We illustrate our scheme by presenting decompositions for various nonlinear 
Hamiltonians including quartic Kerr interactions. Finally, we conclude with two potential exper- 
iments utilizing ofHine-prepared optical cubic states and homodyne detections, in which quantum 
information is processed optically or in an atomic memory using quadratic light-atom interactions. 
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Introduction- Since the proposal of quantum compu- 
tation as a generalization of computer science, an impor- 
tant theoretical challenge has been how to decompose an 
arbitrary gate into a universal set. The corresponding 
theory of discrete-variable decompositions is very exten- 
sive and mostly employs matrix representations of logic 
gates utilizing matrix decomposition techniques In 
contrast to discrete-variable theory, there is not an es- 
tablished method to decompose an arbitrary operator in 
the continuous- variable (CV) regime except the proof-of- 
principle results on universal gate sets in Refs. [ill- In 
particular, Ref. makes use of an exponential operator 
approximation and proves that by employing certain el- 
ementary gate sets (discussed below) one can derive any 
operator up to a certain error. However, none of these 
works intend to present a constructive and efficient de- 
composition recipe. 

The problem of decomposition is intrinsically related to 
the concept of universality. Universality means to have a 
set of operators that allows you to simulate any operator 
on a certain Hilbert space through concatenations of the 
elements of the universal set. For our purpose, achieving 
universality is then equivalent to decomposing, at least 
approximately, an arbitrary unitary exponential operator 
to a set of elementary unitary exponential operators: 

^itH(a,a^ ) _ |giti/fi(a,a*) ^it^H-zia^a^ ) gitjv-H"jv (a,a* ) | 

Here, a and are annihilation and creation operators, 
respectively, and {-ffn} are fixed Hermitian functions of 
mode operators. The coefficients ti,t2, ■■■ are interaction 
times of the Hamiltonians and are functions of t. Thus, 
different concatenations of elements of this set for vary- 
ing interaction times should enable one to simulate an 
arbitrary operator. We assume that we have access to 
arbitrary interaction times for the initial set. 

In our setting, there are two important criteria for CV 
gate decompositions: how systematic and how efficient 
the decompositions are. Here we shall derive methods 



according to these criteria and present a systematic and 
efficient framework for decomposing any given unitary 
operator that acts on bosonic modes into a universal set 
of elementary CV gates. Our general method consists 
of first expressing operators in terms of linear combina- 
tions of commutation operators and then realizing each 
commutation operator and their combinations through 
approximations. We will discuss the efficiency of the 
decompositions and present guidelines to obtain an ar- 
bitrary order of error. For this purpose, we employ a 
powerful technique for obtaining efficient approximations 
(see supplemental material for details). Throughout, we 
use the convention h = 1/2, i.e., the fundamental com- 
mutation relation is [X,P] = i/2 with X = {a^ + a)/2 
and P = i(at - a)/2. 

General Gaussian decompositions- For Gaussian op- 
erators, i.e., second-order operators, exact and finite de- 
compositions to elementary sets are known (here, order is 
defined as the polynomial order of the mode operators in 
the Hamiltonian of a given operator). For example, the 
Bloch-Messiah decomposition allows to decompose any 
second-order operator, i.e., any unitary Gaussian oper- 
ation, to passive linear multi-mode optics, single-mode 
squeezing, and displacement operations 

In Ref. 3 the following set is presented as a single- 
mode Gaussian universal set: {e*^^"^ \ 6**1"^, e**^"^ }, 
and in Ref. similar to the Bloch-Messiah decompo- 
sition, a recipe is given to decompose any single-mode 
transformation of second order to this set with no more 
than four steps. These exact decompositions emerge from 
the fact that mode transformations through second-order 
unitary operations are linear, and thus, one can utilize 
matrix representations and matrix decomposition tech- 
niques. In fact, for Gaussian operator decompositions, 
one can find infinitely many elementary sets and decom- 
positions, and instead of those sets above, one may choose 
the one that suits best the given situation and purpose. 

Universal decompositions- In order to decompose an 
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arbitrary single-mode operator in CV systems, it has 
been shown that, in principle, adding a nonlinear ele- 
ment (of order three or more) to the toolbox is sufficient 
0. In the present work we use the following set: 



{e 



if (X^+p2) 



e 1 ,e " , 



e j. 



(1) 



This set is not unique, and one may use different Gaus- 
sian elements as explained above and a different non- 
linear element. However, this particular set turns out 
to be useful for describing CV quantum computation 
in the one-way model using CV cluster states [1, 01 ■ 
Note that one can simplify this elementary set further by 
omitting the second-order Hamiltonian, since e** ^ — 

27 e 3 e a e e 3 , and usmg the Fourier 

transformation whose action is given in Eq. ([2]). Even 
though this simplification has value from an academi- 
cal point of view, as it reduces the minimal number of 
elementary gates, we are basically motivated by decom- 
posing an arbitrary gate to a set of experimentally acces- 
sible gates. All second-order gates are relatively easy to 
implement and replacing them by third-order Hamilto- 
nians will increase the complexity of the gate sequence. 
Therefore we shall use the overcomplete set ([T]) in our 
decompositions without loss of generality. One may also 
prefer a further extended set depending on a certain ex- 
perimental situation in order to reduce the complexity of 
the decompositions. 

In addition, one can obtain some nonlinear opera- 
tions through unitary conjugation: [/e'*^(°'° _ > 
gitH(;7a(7+,(;7a(7+)+) important unitary conjugation is 
the Fourier transform: 



-if (X^+P^) 



itP" 



(2) 



Employing unitary conjugation, with the set ([IJ, one can 
now generate certain nonlinear gates exactly. For exam- 
ple, e^*^'e'*P'e-**^' = e'*(^-*i^')' , which is a fourth- 
order operator. However, there is only a limited number 
of such decomposable nonlinear operators, and therefore, 
for generality, we will make use of the idea of operator 
approximations . 

Besides the abstract notion of universality [1] how can 
a given unitary exponential operator be decomposed to 
the elementary set ((IJ in a systematic and efficient way? 
Let us first introduce the available toolbox employed to 
achieve a decomposition as efficient as possible (while ef- 
ficiency will be defined and discussed later) . The tools we 
use for CV gate decompositions include Gaussian oper- 
ator decompositions, unitary conjugation, and exponen- 
tial operator relations as well as approximations. Before 
proceeding to the general case, let us demonstrate how 
to realize a particular nonlinear exponential operator us- 
ing the above tools and the set ([T]). A very important 
example is the Kerr interaction operator. It allows to 
convert a coherent state into a cat state M and to realize 



a controlled quantum gate for qubits The one-mode 
self-Kerr interaction, up to a Gaussian element, corre- 



sponds to e'*^'^ 



^rt{x^+x^p^+p'x^+p^)_ In order 



to decompose the Kerr operation to the set ([T]), we would 
first write the full Kerr Hamiltonian as a linear combi- 
nation of commutators and then realize each through 
operator approximations. The following relations, to- 
gether with the Fourier transform ([2]), are enough to re- 
alize this gate up to a phase: = -^[X^,[X^,P^]], 
X2p2 _,_ p2^2 ^ _|[x3^p3]^ rpj^^g^ jg sufficient to 

realize the above commutators and their linear combina- 
tions (see supplemental material for more details). 

Let us now present the general scheme for an arbitrary 
Hamiltonian. Obviously, any single-mode Hamiltonian, 
as a polynomial of bosonic mode operators, consists of op- 
erators of the form cX"^P"+c*P"X"^. We can show that 
any such operator can be written as a linear combination 
of commutation operators. First note that cX"^P" + 
c*P"X"' = Re(c)(X™P" + P"X"^) + iIm(c)[X™,P"], 
and then one can derive the following two identities (see 
supplemental material for a derivation), 

^ -[X"^-\[X\P\ (3) 



X' 



3(m 



j^m pn _j_ p''^ X^ 
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n+l 



rn+l pn-\-l 



(4) 



k=l 



Equation (|3]) is necessary to obtain arbitrary powers of 
X and P operators with the Fourier conjugation ([2]), and 
Eq. (jll) basically prescribes how to systematically decom- 
pose an elementary Hamiltonian to commutation opera- 
tions of orders of X and P and their combinations where 
we can use the tools we have. 

For multi-mode operators one needs an extended ele- 
mentary set including an entangling operation [2| , for ex- 
ample in the optical context it can be the beam-splitter 
operation. Using Gaussian decomposition methods, for 
simplicity, we may assume that we have access to the fol- 
lowing gate without loss of generality, Cz = e^**'^!'^^^. 
For multi-mode Hamiltonians we can again use the sim- 
plifications for a single mode and Eq. (g]) , because of the 
fact that the operators on one mode commute with the 
operators on the other mode. However, for this purpose, 
we initially need to realize the two-mode operations with 
arbitrary powers of X and P in both modes, similar to 
the single- mode relation ([3]). The following relation to- 
gether with Fourier conjugation ^ and Eq. ([3]), is suf- 
ficient to realize the two-mode operations with arbitrary 
powers of X and P, 

A" ®Pi = - [P2^\ [pr\ ^1 ® ^2]] . (5) 

Then, we can use Eq. Q again with the single-mode op- 
erations to realize an arbitrary two-mode expression. As 
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an example, consider the cross-Kerr Hamiltonian (up to 
a Gaussian transformation): (X^ + (g) (X^ + P^)2 = 
XI (g) X| + (g) P| + P2 X| + P2 p2 , In this case, the 
nested commutator [P^j "X" -''^2]] will suffice. An- 

other, even simpler example for decomposing a nonlinear 
two-mode evolution, namely cubic parametric down con- 
version with a quantized pump, is discussed in the sup- 
plemental material. In form of a dispersive interaction, 
it may also be used to mimic a Kerr gate [lo| . 

From an academic perspective, Eqs. ([3]), ([4]), and ^ 
are universal not only for any Hamiltonian, but also 
for any initial universal set with a nonlinear gate dif- 
ferent from X^ because of the well-known equations: 
1^ = -2i[X, P], If = 2i[P, P], where P is a function of 
operators X and P. Thus, any initial nonlinear Hamil- 
tonian can be reduced to a form X™. From a practical 
perspective, however, it is unwise to use Eqs. ([3]) and 
(|4]) for arbitrary initial sets because of a typical increase 
of complexity. Instead, one should derive an optimized 
expression (in terms of the decomposition efficiency, see 
below) utilizing the available tools for every other Hamil- 
tonian and every other universal set. 

Efficiency- Besides having a systematic method, we 
also require the decompositions to be relatively efficient. 
We define efficiency as the number of operators needed 
to realize, in an approximate fashion, a given operator 
with a certain negligible error (note that this definition 
slightly differs from previous ones ^ where efficiency is 
the scaling of the number of operators with respect to 
the error). Let us give a few more definitions. If in 



(6) 



the Taylor expansion of both sides matches for the orders 
of t up to t™, it is called mth-order decomposition [3^ . 
For example, an important case is when C = A + B for 
which we will use the term splitting. Another important 
case is when C — [A, B] which, from now on, we call 
commutation operator. For example, the identity below 
is a well-known second-order approximation for a com- 
mutation operator. It has been used already in quantum 
control flU , discrete- variable quantum computation [l^] , 
CV quantum computation or, in general, Hamilto- 
nian simulation theory [l3| . 



+ f{t^A,B) 



(7) 



It basically says that, for t < 1, the corresponding oper- 
ator concatenation is the same as applying the commu- 
tation operator of A and P, up to some error where the 
dominant term is of the order t^. Now in order to obtain 
more reliable gates, a straightforward and common way 
to improve accuracy is using smaller interaction times, 
t — > t/n, and applying the decomposition times to 
obtain the same interaction time as before. 



Let us call this approach rescaling ([ij] and Refs. 
therein). Besides improving accuracy, rescaling is ab- 
solutely necessary to realize nested commutations. For 
example, one may replace it A by t^[P,A] in Eq. ^ to 
simulate the nested commutation operator e'* [^■[^■-^H — 
^^tB^t^[B,A\^-^tB^-t■'[B,A]J^f'{t^^A,B), and similarly for 
further nested commutations. However, in this identity, 
an approximation of [P,^] is still needed and eventu- 
ally, using again Eq. (O, we obtain an operator whose 
interaction time is of the same order as the dominant er- 
ror term. In order to obtain a reasonable decomposition, 
the order of the dominant error should be smaller than 

. Thus, again rescaling is needed, requiring relatively 
many operators to enhance accuracy. For instance, us- 
ing the approximations described so far, the number of 
operators to approximate a single commutation operator 
with coefficient 0.1 and dominant error term ~ 10"'^ re- 
quires 4000 operators, while for the nested commutation 
operator, with the same values, the number of opera- 
tors will be ~ 10^° (see supplemental material). Hence, 
better approximations are crucial to achieve reasonable 
decompositions. 

What we propose to use is a general and powerful 
method for obtaining higher-order approximations in a 
direct fashion (details on this rather technical part of our 
proposal can be found in the supplemental material; see 
also, for instance, page 393 of Ref. [l^]). In this approach, 
concatenations converge much faster to an arbitrary set 
of commutations and linear combinations of these, reduc- 
ing the number of operators from the order of 10^" to the 
order of tens. 

Experiment- An immediate consequence of our work 
is that it bridges the originally huge gap between ab- 
stract notions of CV decomposition theory and possi- 
ble experimental implementations, for instance, within 
quantum optics. Highly nonlinear quantum gates (such 
as quartic Kerr gates) may then be realized in a de- 
terministic fashion by concatenating tens of quadratic 
and cubic gates. Figure 1 illustrates two experimental 
schemes. They both combine ideas for conditional opti- 
cal state preparation, in order to probabilistically prepare 
and distill high-fidelity cubic phase states ~ Jdxe^*^ \x) 
IB - isj l , with those for temporally encoded (Tol [20| and 



J^[A,B] _ 



fully homodyne-detection-based [2l| CV cluster compu- 
tation. The complication of realizing nonlinear gates is 
shifted offline into the preparation of the cubic ancillae. 

In Fig. la), the conditional state preparer (CSP) is 
linked with a quantum memory (QM) through a classi- 
cal channel in order to signal whenever an offiine state 
is available such that an optical pulse carrying the lat- 
est quantum information is released from the quantum 
memory. After the first application of the optical Cz 
gate between a first input pulse coming from the left 
and a first ancilla pulse coming from the top (emerging 
from the CSP), the C^-transformed input pulse is mea- 
sured through homodyne detection (HD), while the Cz- 
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a) 




b) 



CSP 



spatial encoding using 5.5dB-squeezed light sources [27 1. 

In summary, we proposed a systematic framework to 
decompose arbitrary CV unitaries (from arbitrary multi- 
mode Hamiltonians) into an experimentally realizable 
elementary gate set. Different from previous proof-of- 
principle demonstrations, our treatment brings the ab- 
stract notions of decomposition theory for CV quantum 
computation close to experimental implementations. 

We acknowledge support from the Emmy Noether Pro- 
gram of the DFG in Germany and thank Nick Menicucci 
and Akira Furusawa for useful discussions. 



FIG. 1: Concatenating elementary gates using offline condi- 
tional state preparation (CSP) of squeezed and cubic states, a 
quantum memory (QM), an optical Cz gate, g2»*-^i'?'^2 ^ 
homodyne detection (HD); solid lines represent the optical 
paths. In a), an optical quantum state is released from the 
QM only when the CSP succeeds (dotted line is classical feed- 
forward). In b), quantum information storage and processing 
go hand in hand by employing suitable quadratic light-matter 
interactions between optical pulses and an atomic ensemble. 



transformed ancilla pulse is sent to the quantum memory 
with its quantum state transferred onto the memory. Af- 
ter preparation of the next ancilla pulse, a new input 
pulse emerges from the memory for a second application 
of the Cz gate and so on. This is similar to the scheme 
of Ref. [IQi] for a single quantum wire (corresponding to 
a linear CV cluster chain), but this time including cubic 
states and quantum memories, and excluding non-linear 
measurements such as photon counting (except for the 
CSP). Note that storage of non-classical states has been 
already demonstrated experimentally [l^, [l^ . 

The scheme in Fig. lb) uses an atomic ensemble and a 
quadratic light- matter interaction [23|. The light- matter 
interaction is always delayed until an appropriate offline 
state is available. Quantum information can be stored 
and processed at the same time within the ensemble [2^ 
by swapping the optical and atomic states after every 
interactiori, simply using additional quadratic interac- 
tions m, [i^. As opposed to Refs. [ll, [l^, the opti- 
cal state is either a squeezed vacuum or a conditionally 
prepared cubic state, inserted into the temporal cluster 
chain whenever needed [2lj. Thus, we make explicit use 
of the atomic memory to compensate for the probabilistic 
optical ancilla-state preparations. Only with these cubic 
ancillae is it possible to perform universal gates solely 
through homodyne detection 2l|. In principle, our de- 
composition method would then determine what homo- 
dyne local-oscillator angles to choose in order to realize a 
given Hamiltonian in a deterministic fashion. Currently 
available squeezing levels of almost 13dB are, in principle, 
sufficient to perform up to 73 elementary teleportations 
in a nonclassical fashion; four such elementary Gaussian 
gates have been implemented already in a fully optical. 
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[32] The approximation quality does not only depend on the ing operators. In the present work, the focus is on the 

decomposition order, but also on the actual values of the decomposition order, omitting the absolute errors, 

gate interaction times and the norms of the correspond- 

SUPPLEMENTARY MATERIAL 
DERIVATION OF EQUATION (4) 





pnj^m 




j^mpn 


_ pn—l^j^mp 




j^m p n 


_j_ pn—lj^m- 


-1 _ pn-2 pj^m p 


j^mpn 




-1 pn — 2j^mp2 i pn—2j^m—lp 

2 


j^mpn 




— 1 _|_ pn—2 p pn—2j^'mp2 



n-l 

_ im \ - pk ym-l pn-k-l 
fc=0 



n-l 

im 



[X"'' P"] = ^ ^ pk -^m—l pn—k—l _|_ pn—k—lj^m—lpk 

^ fe=0 

. n-l 

^ ^ ^(J^"^~"^P^ [^X^"^""^ jDfcj^jDn— fc— 1 _|_ pn—k—lQ^^m—1 pk"^ _|_ /^'''^X'"^"-^^^ 

^ fe=0 

^ ^ ^j^m— l^n— 1 [J^^""*^ pk'^pn—k—1 _|_ j^n— fe— 1 j^j^m— 1 J?''^! _|_ pn—lj^m — 1^ 

^ fe=0 
n-l 

^ ^ ^j^m— l^n— 1 _|_ pn—lj^m—1 _|_ ^j^n- fe — 1 [_^^~1 -^^]]) 

^ k=0 

/n-2 \ 
_ imn ^j^m-lpn-l _|_ pn-lj^TO-l'J _|_ ^ j "y^Jpn-k-l JJ^"»-1 p^]] j 



\k=l 



As a result: 



4? 1 

x™p"+p"x™ = -7 — - — -[x'"+\p"+M y[p"-'=,[x™,p'=ii 



(n+l)(m+l)' , - 

Note that for n = 1, the summation term in the identity above is zero. Also, due to the Jacobi identity, we have 
l^pn-k^ [X™, P'^]] = [P'', [X™, P"^*^]], and this may also lead to some simplification depending on the value of n. 

NECESSITY OF GOOD APPROXIMATIONS 



Here, we shall illustrate the necessity for having better (than any commonly used) approximations. The approxi- 
mations that we employ are explicitly introduced below in this supplemental material. 
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For instance, for a nested commutation approximation of an interaction strength 0.1 and a dominant error term 
of 10~^, we would need 73 operations corresponding to an eighth-order approximation. For comparison, note that it 
is also possible to use Lloyd's method [13]. Lfoyd's idea was originally intended as a proof of principle, but it has 
also been used in the literature as an approximation tool. The required number of operators using Lloyd's method is 
evaluated below. The notation N[x] is used to indicate that a number of x operators is required. For the commutation 
operator, we have 

+ (8) 

^t[A,B] = ^[4 X n^l + / f — 
V n 

while for the nested commutation, we obtain 



^ ^^B±^A±^-^B±^-^Ai 



irB^inlB.A]^-^i-B^~^lBA] 



i^e'B^^^A^^-^B^^-^A^\^ ^~.^B ^^-W ^ ^-^A^ ^^B ^ ^^A^\^ \ _^ f' (L) ^ f 



^jtlB,lB,A]] ^ ^i2^ ^3] _^f:f _^fft 



Thus, using these approximations, the number of (elementary) operations will be of the order of 10"'^° for an interaction 
strength of 0.1 and a dominant error term of 10^^. 

EXAMPLE: ONE-MODE SELF-KERR INTERACTION GATE 

For the decomposition of a Kerr gate with coefficient 0.1 and dominant error term 10"'^, we need to approximate the 
following operators: d = -|[^^ [X^,P^]], O3 = ~f[X^,P^]. The Kerr gate is then : 6^01(01+02+03). Heje Oa^is 



the Fourier transform of Oi, i.e., e'*02 — Fe**Oi^t. First, using second-order three-party splitting (see Refs. 14l.l28|). 
we split into separate elements: 

^^0.1iO^+O,+Os) ^ ^^0.05O^ ^^0.05O, ^^0.1Os ^^0.05O, ^^0.05O^ ^ ^ j [0^02, O3) + . . . (10) 

and then we insert the approximations for the commutation and the nested commutation operators (see below). 
Approximations for these operators can be calculated considering the necessary approximation order. For the nested 
commutation operator with coefficient 0.05 x 2/9 and dominant error term smaller than lO"'^, we need the fourth-order 
approximation which corresponds to 9 operators, 

g-i0.05 X I [X^ , [X^,P^]] _ g- iO. 11157p2 gj0.02231X^ ^ia.Q2231P^ g- j0.02231X^ g-i0.0223ip2 

g-»0.02231Jf^g»0.02231p2g.0.02231X^g-^0.11157p2 55326 x lO"^ X f{X',P^) + ... 

Similarly, the commutation operator with coefficient 0.1 x 4/9 and dominant error term smaller than 10^'^ also requires 
the fourth-order approximation, this time corresponding to 10 operators, 

-0.1 X I [X^ ,P^] _ -i0.25298X^ i0.210819P^ i0.01918X^ -i0.28476P^ i0.36163Js:^ 

gj0.36053P^gi0.12861X^g-i0.05805P^g-i0.25644X^g-i0.22853P^ _|_ q 4J^643 X 10"'^ X f{X^ P^) + 

As we have to apply Oi four times and O3 only once, we require 4 x 9 -I- 10 = 46 elementary operations. Every single 
Oi operation as well as the O3 operation each consume 10 additional Fourier transforms in order to switch from 
the elementary X gates to the necessary P gates [here, ~ 10 indicates that some Fourier gates cancel in the sequence 
PH)) when switching between Oi and 02]- As a result, summing up, we will need 9 -I- 10 -I- 9 -I- 11 -I- 9 = 48 extra 
Fourier gates, in addition to the 46 elementary quadratic and cubic X gates. In total, this leads to 94 elementary 
operations from the universal gate set needed to simulate a Kerr interaction of size 0.1 with errors scaling smaller 
than 10^"^. Without using these more powerful approximation techniques (as explained in detail below) and using 
instead standard techniques such as the well-known Trotter formula for splitting and the standard approximations for 
the commutation operator ([5]) and the nested commutation operator ([5]), we would need '^10^ operations to simulate 
this Kerr interaction gate with the same precision. 
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EXAMPLE: PARAMETRIC-DOWN-CONVERSION HAMILTONIAN WITH QUANTIZED PUMP 

As a second example, we discuss the decomposition of the Hamihonian evolution of parametric down conversion with 
a quantized pump. This example illustrates that one can sometimes find exact operator relations and thus reduce the 
use of operator approximations to some extent. The operator to be decomposed is: f^i-ii^i^^-Pi^-2+XiPiP2+PiXiP-2) ^ 
where the subscripts denote the two optical modes. As described, we first divide the Hamiltonian into its simplest 
Hermitian parts {^e^^^-^^ , e'^i , e'(^i^i -^2+^1-^1^2)1 through splitting approximations. Now the operators e^^i^^ and 
giPi X2 realized through exact operator relations instead of operator approximations: 

Here, k and a are some free parameters and e^^^^^ ^e^^^-^^ are equivalent to e'^i^^ up to Fourier transforms. The 
only remaining operator, g'^^i^i^^+^i^^i^^), can be realized through a commutation approximation of [Pl^XiP2\- 



APPROXIMATING COMMUTATION AND NESTED COMMUTATION OPERATORS 



Introduction 



In this section, we shall explain a very powerful method that we propose to employ for deriving the necessary 
approximations in our scheme. We also present explicitly some efficient approximations which are used throughout 
the paper. 

Let us first define the problem. Suppose we have access to a finite set of non-commuting exponential operators, 
jgt gt B ^ gt c ^ _ 1^ where A, i?, (7,... are arbitrary operators. Through different and finite concatenations of these 
operators, our aim is to obtain another operator, {e**^}, approximately, i.e., up to a certain error, 

M 

e*o= J|e*'^e*'^e*"'5... (12) 

i 

where t,ti,t[, ... are scalars. Also, we further assume that ... are not fixed, i.e., they are adjustable parameters 
being functions of t as a result of the approximation problem. If the Taylor expansion of both sides of Eq. 
matches for the orders of t up to t™, then the right hand side of the equation is called rn'th order approximation for 
the operator e**^, and the problem is to find a method for this approximation, i.e., to find the appropriate values for 

^ f' 

i/j, tj, — 

In the present context, we are interested in the case where O is an arbitrary nested commutation operator comprised 
of A and B: [., [., [., ...[A, -B]]]], which is to be approximated through the set {e*"^, e* ^}. Throughout this section we 
call an exponential operator with exponent [A, B\ commutation operator, and similarly for the nested commutators. 
In general, we omit hats to indicate operators except when confusions are possible. 



Employed method 

The method we employ combines various elements of existing splitting and approximation techniques (see, for 



instance, page 393 of the splitting review Ref. (15| ) . Our execution of these methods relies upon the generalized Baker- 
Campbell-Hausdorff series (gBCH) which is basically the logarithm of an arbitrary exponential operator concatenation 
involving more than two operators; an extension of the standard BCH, as the name indicates. It is possible to 
calculate the gBCH in various ways, for instance, by repeatedly utilizing formulas for the standard BCH. With the 
aid of computer software, this is straightforward. In the present work, however, we shall employ a very convenient 
calculation technique proposed by Reinsch [29.]. 

For any necessary orders of the operator approximation, we use two calculation steps. In the first step, we employ a 
brute force and hence efficient method for approximations. In the second step, we obtain higher-order approximations 
using the results of the first step. This second step will be necessary in order to be able to achieve orders of 
approximations otherwise unattainable through the first step alone. 

In order to derive a certain order of approximation, we start with a specific operator concatenation like the one in 
Eq. We calculate the logarithm of this concatenation through the gBCH and this results in a linear combination of 



8 



a set of operators. The coefBcients of these operators correspond to polynomials of the input parameters. We then solve 
all the polynomial equations in order to eliminate the undesired operators in the linear combination. The solutions 
of the polynomial equation set give us the proper values of the input parameters for a particular approximation. 
Hence, the whole problem of obtaining a specific order of approximation is reduced to a computationally well defined 
problem, namely that of solving polynomial equations. 

For example, consider a fourth-order approximation of the commutation operator e* ["^-^l through the operators 
jgit gzt Here, i is the square root of -1. We start from a concatenation of alternating e** and e** ^ operators. 
The concatenation consists of ten operators, for which the reason will become apparent below. Thus, we have 

log(e-i"'^e^^"^e^^^*^...e=i«^*^) =htA + f2tB + f3t^[A, B] + /4t'[^[A, B]] + hv'lB, [A, B]] 

+ ht^[A,[AM,B]]]+ ht^[BMA.m]+ ht^[B,[B,[A,B]]] + ... 

where the linear dependent operators such as [B,^] are omitted. The eight coefficients, {/i, /§}, are polynomials 
of the input parameters, {ci, cio}. To these parameters, we assign the values ci = 1.2 and C2 — —1 to ensure a 
real solution set. Then we solve these polynomials simultaneously with the conditions — 1 and Jk = 0, Vfc 7^ 3. As 
a result, we obtain sets of appropriate values for the input parameters. Consequently, the left hand side of Eq. (fT3|) 
can be considered a fourth-order approximation for the commutation operator. More specifically, for the first step of 
the calculation, we need to take into account the following points. 

Constructing real solutions- the number of operators in the concatenation should always match the number of 
polynomials to be solved. However, the solutions might be all complex, which is useless for physical scenarios. In this 
case, one has to add more operators to the concatenation to be able to adjust the input parameters and obtain an all 
real-valued solution set. This is why in the above example, we added two more operators to the concatenation and 
chose the specific values 1.2 and -1. We choose the number of extra operators and the assigned values by trial and 
error. 

Removing linear dependency- one should reduce the polynomials for simplification. The polynomial orders depend 
on the order of the operators that belong to the polynomials (for instance, in the above example, /a is a second- 
order polynomial, but and /a are third-order polynomials), and for every order there exist many linear dependent 
polynomials. Especially for higher orders, these will consist of thousands of elements. Thus, in order to reduce the 
complexity, one should eliminate the linear dependent polynomials. For example, a concatenation series comprised 
of 14 operators, for a fifth-order approximation, in general, gives rise to 51 polynomials, but, after the elimination of 
the linear dependent ones, this reduces to 14. 

Solving polynomial equations- as an algorithm for polynomial equation solving, we do not recommend to employ 
a conventional analytical method such as those standard methods based upon Grobner bases or other standard 
commercial computer software. These approaches hardly help for obtaining solutions of orders higher than four. 
Instead, for polynomial equation solving, we use the so-called homotopy continuation method fsO] for the first step of 
the calculation scheme (as well as for the second step, as described later). 

The homotopy continuation method is an iterative numerical method. Even though it has its own complications 
and subtleties, we shall only use the basic principle of that method. We employ it for iterating Newton's polynomial 
solving algorithm. Newton's algorithm works only provided one can guess the roots of polynomials very close to their 
actual values. Homotopy continuation is basically a general idea to circumvent this kind of problems by just starting 
from a simple toy problem and slowly converting it into the real problem while tracing the solutions. In our case, 
we start with a set of polynomial equations the roots of which are obvious, and make small changes on these initial 
polynomials. We then solve the new equations using the results from the previous polynomial solving step as the 
starting point for Newton's method. Further, we iterate this process such that at some point we get the original 
polynomials that we wanted to solve, and in addition some of their solutions. Besides the overall time efficiency of 
this approach, the method allows us to deal with each solution separately, thus, eliminating the need of finding all 
the solutions. This property indirectly improves the efficiency dramatically. 

We can obtain solutions for up to fifth-order approximations and with a little more expertise on polynomial equation 
solving and patience, it should be possible to obtain even sixth-order approximations. However, in order to obtain 
higher-order approximations, we believe that the first step alone is not sufficient. Therefore, we propose to employ 
yet another step in order to be able to go to higher orders 

In order to circumvent the problem of having an intractable calculation, we have to employ another step. In this 
second step, to derive higher order approximations, we now use a concatenation of those approximations that we 
obtain through the first step. This second step is similar to the Suzuki- Yoshida method [SUM [HI- The difference 
is, instead of improving step by step and obtaining every higher order recursively, here, by utilizing the gBCH, we 
obtain a by many orders higher approximation in a single step. 
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Let us turn again to our example of commutation operators. Suppose we want to obtain a ninth-order approximation 
for a given commutation operator. In this case, we can obtain the fifth order with the first step, corresponding to 
a approximation for an operator of the form e* ^g-j- denote this as Q^it). We then take a 

concatenation series of two such operators Q5{t) and Q^^{t). Obviously, the approximation of Q^^{t) is the reverse 
ordered approximation of Q^i—t), and hence is an operator of the form g-t^[A-B]-t'^F-t''G-... ^ rpj^g concatenation we 
use is Y[iQ5{Eji)Q5^iPi^) where pi and p'^ are parameters that we have to find (this kind of concatenation is also 



used in Ref. 31| for a recursive approximation in the splitting problem). We again employ the gBCH to obtain the 
resultant operator combination and eliminate the operators from order six to nine. For this particular example, we 
need to solve seven polynomials, but due to the constraint of a real solution set, we use nine operators. When we 
rewrite each fifth-order operator in terms of their approximated form, in total, we need 144 operators. Obviously, this 
method also heavily depends on the efficiency of polynomial solving. 



Examples 

In this section, we present some commutation and nested commutation approximations. Although for each case 
there exist many solution sets, we always present only one of these and we only give the explicit numbers up to six 
digits out of fifteen digits. Also, note that here some parameters are indicated as zero in order to avoid confusion in 
the general concatenation ansatz that we used. 



Commutation operator 
Here, the operator to be approximated is e* ["^-^l and the general ansatz we use is: 



M 

i 

A fourth-order and fifth-order approximation is presented in table [II and table [III respectively. 



i 


Ci 


/ 

Ci 


1 


1.2 


-1 


2 


-0.090992 


1.350762 


3 


-1.715364 


-1.710162 


4 


-0.610065 


0.275377 


5 


1.216422 


1.084021 



TABLE I: fourth-order commutation operator 



i 


Ci 


Ci 


1 


1.2 


-1 


2 


-2.121328 


-1.680681 


3 


-0.0199879 


1.602050 


4 


2.025319 


2.743677 


5 


0.061217 


-1.170812 


6 


-3.235890 


0.018735 


7 


1.519023 


-0.989125 


8 


0.571646 


0.476156 



TABLE II: fifth-order commutation operator 
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i 


Pi 




1 


0.8 


-1 


2 


1.427601 


1.727306 


3 


1.816745 


1.700274 


4 


1.371191 


-1.693846 


5 


-1.698488 






TABLE III: ninth-order commutation operator approximation through fifth-order approximation 

We can then employ the second step of the calculation scheme to generate a ninth-order approximation through 
fifth-order approximations, as shown in table Hill For this purpose, the concatenation ansatz we use is: 



l[Q5{p^t)Q5\p^t) 



(15) 



Nested commutation operator 

Here, the operator to be approximated is e**^['^'["^^-^ll and the general ansatz we use is that of Eq. (HH). Respectively, 
in table IIVI and in table |Vl fourth-order and fifth-order approximations are presented. 



i 


Ci 


1 


1 





0.5 


2 


-1 


-1 


3 


1 


1 


4 


1 


-1 


5 


-1 


0.5 



TABLE IV: fourth-order nested commutation operator 



i 


Ci 


/ 

Ci 


1 





0.5 


2 


-0.912433 


-1.000141 


3 


2.439891 


-0.531940 


4 


0.184788 


2.000998 


5 


-0.477395 


-1.659152 


6 


-0.761744 


1.465192 


7 


1.181123 


-1.312200 


8 


-1.654230 


0.537205 



TABLE V: fifth-order nested commutation operator 



In table IVTl a ninth-order approximation is presented using the ansatz of Eq. (jl5p . 

Given an initial approximation, it is also possible to use the symmetry argument of Suzuki IJ, |31| and construct 
recursively symmetrical approximations for the nested commutation operator, since e** [^'[^'^lle'(~') [^.[-4,5]] _ j 
Hence, a sixth-order approximation for nested commutation can be constructed through fourth orders using 3x9 — 2 = 
25 operations. 



i 


Pi 




1 


1 


1.2 


2 


-0.948048 


-0.830428 


3 


1.049237 


-0.887661 


4 


1.358200 


1.329989 



TABLE VI: ninth-order nested commutation operator through fifth-order approximation 



